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We show that spherical truncations of the 1/r interactions in models for water and acetonitrile yield very 
accurate results in bulk simulations for all site-site pair correlation functions as well as dipole-dipole correlation 
functions. This good performance in bulk simulations contrasts with the generally poor results found with 
the use of such truncations in nonuniform molecular systems. We argue that Local Molecular Field (LMF) 
theory provides a general theoretical framework that gives the necessary corrections to simple truncations in 
most nonuniform environments and explains the accuracy of spherical truncations in uniform environments by 
showing that these corrections are very small. LMF theory is derived from the exact Yvon-Bom-Green (YBG) 
hierarchy by making physically-motivated and well-founded approximations. New and technically interesting 
derivations of both the YBG hierarchy and LMF theory for a variety of site-site molecular models are presented 
in appendices. The main paper focuses on understanding the accuracy of these spherical truncations in uniform 
systems both phenomenologically and quantitatively using LMF theory. 
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I. INTRODUCTION 

Spherical truncations of the Coulomb interactions present in typical molecular models such as CHARMM [Ql |2l and AM- 
BER have long been used to keep computational cost in check. This cost in the simulation of large biomolecules is com- 
pounded by the use of explicit water models containing point charges to describe the hydrogen-bonding network and dielectric 
behavior of the solvating water molecules. Since traditional particle-mesh Ewald sum treatments of Coulomb interactions do 
not scale well in massively parallel simulations, a computationally compelling case can be made for the use of spherical trunca- 
tions |4j]. However, spherical truncations have been shown to be clearly wrong when applied naively in a variety of nonuniform 
environments iHI^. For this reason, the use of short-ranged truncations of 1/r interactions is typically viewed as an unjustified 
approximation. 

There have been many attempts to place the use of spherical truncations of 1/r on a more solid theoretical footing, including 
site-site reaction field methods 101, Wolf summation flS, and isotropic periodic summation lflol - [l2ll . Despite this work, the 
virtues and defects of spherical truncations of 1/r in various applications remains a subject of ongoing debate in the current 
Uterature 

Our approach, local molecular field (LMF) theory lfl6l [TtIi . uses an effective single particle potential to account for the 
averaged effects of the long-ranged interactions neglected in typical spherical truncations. It gives a theoretical basis for the use 
of simple truncations in some cases, and also provides a physically suggestive path for correction when such truncations fail. 
Moreover, recent work has established a very efficient and accurate numerical method to determine the effective field in LMF 
theory using a simple linear response approach ifisll . 

LMF theory for general nonuniform fluids is derived from the exact statistical mechanical Yvon-Born-Green (YBG) hierar- 
chy lfT9ll20ll by making two physically-motivated and well-founded approximations. These rely on the ability of well-chosen 
truncated potentials to yield accurate nearest-neighbor correlations and on the corresponding slowly-varying nature of the re- 
maining long -ranged parts of the full potential iflTll . Previous work has shown that LMF theory corrects two well known failures 
of spherical truncations of 1/r interactions: 



• simulations using LMF theory yield correct charge density profiles for water confined between two walls 112 lH and for ions 
confined between charged plates ll22ll . and 

• simple analytical corrections derived via LMF theory result in accurate energies and pressures for uniform ionic and 
molecular systems lfT6ll23ll treated with spherical truncations. 

In this paper we employ LMF theory to illustrate and explain why spherical truncations of 1/r can often be applied very 
successfully for determining the structure and thermodynamics of uniform molecular systems. When LMF theory is applied to 
charge-charge interactions, all 1 /r interactions are split into short and long ranged parts vo{r) and f i (r), such that 

1 _ crfc(r/cr) crf(r/cr) 

- = Vo(r) + vi(r) = \ . (1) 

r r r 

Here vi{r) is the electrostatic potential from a unit Gaussian charge distribution with width a, and uo(r) corresponds to the 
potential from a point charge surrounded by a neutralizing Gaussian charge distribution lfl6ll . Thus wo(^) vanishes at distances r 
much greater than the "smoothing length" <t and at small distances the force from vo{r) approaches that of the bare point charge, 
so uo(r) can be though of as a "Coulomb core potential". 

In the simple strong coupling approximation (SCA) to the full LMF theory, we assume that all effects from the long-ranged 
interactions due to vi{r) may be neglected. Thus the SCA is in essence a spherical truncation where all 1/r interactions are 
replaced by the short-ranged wo(r), with a setting the scale for the truncation distance. In Section|II] we emphasize the accuracy 
of the SCA for uniform molecular systems, presenting results for SPC/E water and acetonitrile (CH3CN), including the effect 
of varying the range of the short-ranged truncation of 1/r as represented by a in equation These results can be appreciated 
independent of the underlying LMF theory discussed in the rest of this paper. 

Furthermore, we demonstrate that spherical truncations can lead to highly accurate dipole-dipole correlations in uniform 
molecular systems. This surprising result is in sharp contrast to findings by Nezbeda, who used a different molecular-based 
truncation scheme ll24ll . and we shall explain our success later using the full LMF theory. Then in Sections Hill and [TVl we 
formulate LMF theory for bulk uniform site-site molecular fluids and discuss the success of these spherical truncations and 
the neglect of long-ranged interactions using the LMF theory framework for the simpler bulk water system. The form of the 
derived LMF equation and the necessary approximations make clear why spherical truncations can often give accurate structure 
in uniform systems, despite their invalidity in nonuniform systems. 

Detailed derivations of LMF equations for various molecular models are discussed in complementary appendices. Here we 
build on previous derivations of LMF theory for simple atomic fluids, and focus in particular on the derivation of the LMF 
equation for a uniform system of site-site molecules, described by the Hamiltonian 

N ^ N N n n 

U = Y, u:m{B..) + 2 E E(1 - E E - !)• (2) 
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FIG. 1. Schematics depicting tlie geometry and site labels of the water model 12511 and acetonitrile model I26h used in this paper. 



Here Rj describes the positions of all sites within a molecule i connected by a generalized bonding potential ojM(R-i), and 
Ma^drl"'' — rj^"* I) describes the general pair interaction between two sites a and ^ on two different molecules i and j as insured 
by the Sij term in equation (|2]). In AppendixlAl we present a derivation of the LMF equation used in previous work for small site- 
site molecules in a general external field. Then in Appendices IB] and ICl we present the notationally more complex derivations 
of both the exact YBG hierarchy and the LMF equation for a uniform fluid composed of these site-site molecules. Finally, in 
Appendixini we present an abbreviated derivation for larger molecules described by CHARMM- or AMBER-Iike Hamiltonians, 
thus supporting the validity of our conclusions for systems composed of much larger molecules. 



II. STRONG COUPLING APPROXIMATION (SCA) SIMULATIONS OF WATER AND ACETONITRILE 

We present structural results for the simulation of two different small site-site molecular models shown in Fig.[T] 

• SPC/E water jzsll, a rigid molecular model of a hydrogen-bonding fluid, and 

• acetonitrile, an AMBER-Iike flexible molecular model ll26ll of a strongly dipolar fluid. 

These models, along with annotation used for each site, are shown in Fig.[T] 

For the water simulations, we present results for simulations of 1728 SPC/E water molecules in a cubic box of side length 
37.27 A using the DLPOLY2.16 simulation package IztIi . The system of water molecules was equilibrated for 500 ps at 300 K 
using a Berendsen thermostat with a time constant of 0.5 ps and a timestep of 1 fs. Data was collected over the subsequent 1.5 
ns. Cutoff radii for the potential ranged from 9.5 A up to 13.5 A for varying choices of a in ?;o('''). The SCA simulations are 
compared to simulations using Ewald summation with a = 0.3 A~^ and kmax = (10, 10, 10). 

We have previously shown that the SCA with a a of 4.5 A gives a highly accurate O-O pair correlation function for SPC/E 
water ll2lll . In Fig.|2a-c), we show the pair correlation functions for all site-site pairs using the SCA with a ranging from 3.0 A 
to 6.0 A. In all instances, the g{r) are in excellent agreement with results of the full system determined using Ewald sums. In the 
plot of gHH(?'), the curves for each cr choice are displaced by 0.2 in order to emphasize that all plots contain multiple choices of 
vo{r) while yielding essentially the same correlation functions on the scale of the graph. 

These data illustrate the important point that cr is a consistency parameter rather than an empirical fitting parameter lfT6l[l7ll . 
Thus the mean field averaging leading to LMF theory become highly accurate for any choice of a greater than a state dependent 
minimum value (Tmin, typically of order a characteristic nearest neighbor distance. For SPC/E water tTmin is about 3 A, the radius 
of the Lennard- Jones (LJ) core on the water molecule. Any smaller cr would clearly yield a short-ranged system that does not 
accurately describe the oxygen-hydrogen charge correlations on neighboring molecules that compete with the LJ core repulsions 
in forming hydrogen bonds, and indeed poor agreement is found at smaller cr. 

We also carried out molecular dynamics simulations of bulk acetonitrile at two very different states, a high density liquid 
very near liquid-vapor coexistence at 298 K and a lower density system at 550 K. We used a six-site model with flexible 
bonds developed by Nikitin and Lyubartsev |f26ll in which intermolecular potential parameters have been optimized for better 
consistency with experiments. In order to simulate at appropriate bulk densities, an initial configuration of 864 molecules in 
a cubic box is equilibrated for several hundred picoseconds (ps) in the NPT ensemble using the Nose-Hoover thermostat and 
barostat until the volume has equilibrated. The low temperature system has a simulation box length of 42.2 A. The dilute system 
at 550 K has a density one-third that of 298 K and is further equilibrated in the NVT ensemble for several hundred picoseconds. 
The cutoffs of the Lennard- Jones interactions were set to 15 A. When the SCA is employed, the cutoffs for vo{r) were 15 A for 
cr = 4.5 A, 21 A for cr — 6.5 A, and 30 A for a = 8.5 A. When Ewald summation was employed as a benchmark, the cutoff for 
the real space interactions was set to 15 A, and a ~ 0.26 A^^ with kmax = (15, 15, 15). 

Results for acetonitrile site-site pair correlation functions are shown in Figs.[3]and|4] These figures focus on the pair correla- 
tions at both 298 K and 550 K between a nitrogen site (YN) and all four molecular sites on another molecule. The remaining six 
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FIG. 2. g{r) for each site-site pair of SPC/E water, as well as {cos 6) (r). All plots display the correlations determined via both Ewald 
summation (Full) and spherical truncation via LMF theory. The smoothing length a ranges from 3.0 A to 6.0 A in all plots. In the plot of 
(/HH(r), the curves for each a choice are displaced by 0.2 in order to emphasize that all plots shown contain multiple choices of vo{r). Insets 
in (a) and (d) focus on the region near the peak height, where small errors in the a = 3.0 A curves are just visible. 



intermolecular site-site pair correlations are described just as accurately, and are not displayed for brevity. For the high density 
room temperature system, both a shown yield quite accurate results. Note that the <t values are comparable to those used for 
water, despite the greater size of the acetonitrile molecule. Very poor results (not shown) were obtained with use of a too small 
a = 2.5 A as would be expected. 

For the higher temperature, lower density system, a — 4.5 A performs poorly and is not shown, a = 6.5 A is markedly 
improved, and only the largest cr of 8.5 A yields high quality agreement with the Ewald simulation. This is expected from simple 
scaling arguments since the typical nearest neighbor distance is larger; multiplying (Tmin ~ 4.5 A for 298 K by the requisite 
increase in interparticle spacing at lower densities yields 6.5 A. The need for a somewhat larger (Tniin is likely a result of the 
increasing relevance of more extended conformations of these molecules at lower densities and higher temperatures. Both the 
water and acetonitrile results show that spherical truncations are quite good in bulk fluids, given a sufficiently large truncation 
radius. This is phenomenologically well established in the literature. 

The strong agreement of all the acetonitrile site-site correlation functions, given a sufficiently large a, suggests that the 
angular correlations between these molecules are also accurate, for otherwise many of the unusual functional forms would not 
be reproduced with fidelity. Thus we also examine angular correlations for both water and acetonitrile. 

Fig. |2d) shows the excellent agreement of dipole-dipole correlations in SPC/E water simulated using the SCA with those 
correlations in the full Ewald simulations. Here we plot the average cos between water dipoles as a function of separation 
distance r between the centers of mass of two water molecules. Such good agreement is not a consequence of the relatively 
compact nature of the water molecule. Shown in Fig. |5] are plots of {cos 6) (r) for the acetonitrile system at each temperature. 
We again find quite good agreement between the angular correlations in the full Ewald system and our short-ranged systems. 
This agreement follows the trends found for the simpler site-site distribution functions, with a larger cr needed for the low density 
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FIG. 3. Pair correlation functions for the nitrogen site YN on acetonitrile at 298 K as a varies from 4.5 A to 6.5 A. As before, tlie plot of 
gYN-Ycir) displays the g{r) displaced by 0.2 but all plots display results for Ewald summation (Full) and all choices of a. 



higher temperature state. 

We beheve that the excellent agreement of the dipole-dipole correlations in these spherically truncated fluids is a direct conse- 
quence of the general validity of LMF theory and the accuracy of the strong coupling approximation in uniform environments, as 
we describe in the following section. Nezbeda has previously reported poor results for dipole-dipole correlations in short-ranged 
systems where the determined g{r) were accurate ll24ll . The crux of the difficulties with his chosen cutoff scheme was defining 
these cutoffs on a molecular basis, rather than a site basis. This leads to neglected potentials which actually rapidly vary near 
the cutoff radius, counter to one of the important assumptions of LMF theory as discussed in the next section. Takahashi and 
coworkers |[28il studied the effect of cutoff radii in the isotropic periodic sum approach on various properties of water and found 
for (cos 9) (?') that deviations in this property were minimal and equivalent for cutoff radii greater than 16 A. This cutoff radius 
for IPS can be compared to the cutoff radius of 13.5 A used for <t = 6.0 A in this paper We take their observed "saturation" 
in errors beyond a given cutoff radius as an indication that their spherically-truncated potential satisfies the necessary conditions 
for the validity of LMF theory. 



III. LOCAL MOLECULAR FIELD (LMF) THEORY FOR SITE-SITE MOLECULES 

LMF theory for a general nonuniform system prescribes a mapping from the system of interest where all particles interact via 
their full intermolecular potentials in the presence of an external field to a "mimic system" where particles interact via short- 
ranged truncations of their intermolecular potentials but in the presence of a an effective or restructured field. The restructured 
field accounts for the effects of the long-ranged components of the intermolecular interactions using a mean field average. Far 
from being a simplistic mean-field ansatz, LMF theory has been shown to be strongly rooted in statistical mechanical theory, 
and based on physically-motivated approximations that are well-justified for dense fluid systems. 



6 



2 

Syn-ct 
1.8 

1.6 

1.4 

1.2 

1 

0.8 
0.6 
0.4 
0.2 




2 

SYN-YC 
1.8 

1.6 

1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 







. Full 




\ a = 6.5 A 




\ a = 8.5 A 







2 4 6 8 10 12 14 
r(A) 

(a)gYN-CT('') 









\ " _ ^5=8:5 



















2 4 6 8 10 12 14 
r(A) 

(c)gYN-Yc(r) 




2 4 6 8 10 12 14 
r(A) 

(b)gYN-Hc(f) 




2 4 6 8 10 12 14 
r(A) 



FIG. 4. Pair correlation functions for the nitrogen site YN on acetonitrile at 550 K as a varies from 6.5 A to 8.5 A. As before, tlie plot of 
gYN-Ycir) displays the g{r) displaced by 0.2 but all plots display results for Ewald summation (Full) and all choices of a. 



LMF theory for charged systems takes a particularly simple form when charges alone are treated via LMF theory using a single 
a, where the results can be exactly re-expressed in terms of the total charge density and a restructured electrostatic potential. 
Based on the splitting of 1 /r defined in equation ([T]i, each pair potential (r) in equation (|2]i may be decomposed as 

Ua^{r) ^ UQ a^ir) + ^^^Vi{r), (3) 
e 

where UQ^a^{r) contains all the non-electrostatic Lennard- Jones-like pair interactions as well as a uo(r) contribution appropri- 
ately scaled by charge and the dielectric constant e. The crucial feature of these two potentials vo{r) and ui(r) for the validity 
of LMF theory is that tr is chosen so that vq (r) contains all relevant strong Coulomb forces between nearest neighbors and that 
vi (r) is consequently slowly-varying over the range of strongest correlations between those neighbors. A more careful statement 
of the relevant approximations may be found in AppendixlAl 

Previous work focused on nonuniformity such as confining walls, using the Coulomb LMF equation lfl7ll2ll] 

Vn{r) = Vo(r) + i 1 dr' p%,Jr')v, (|r - r'|) , (4) 

where Vq results from the convolution of the fixed charge density with vq (r), and p'j^ includes both the fixed and mobile charge 
densities. Note that Vr and are implicitly functionals of one another, so this is a self-consistent equation. 

Since vi (r) is the electrostatic potential arising from a Gaussian charge density with width cr, equation (HI suggests that the 
restructured external potential Vn may be understood as an electrostatic potential containing the full impact of fixed charges 
and the Gaussian-smoothed impact of mobile charges. In Appendix lAl we present a derivation of the LMF equationj^or site- 
site molecular models as used in previous papers. Equation (01) is identical to that for mixtures of charged species |17], and a 
derivation for small site-site molecular models requires only one further approximation, requiring that intramolecular correlations 
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FIG. 5. Angular correlations represented by (cos 9) (r) for acetonitrile at 298 K (left) and 550 K (right) as a varies. 



are well represented by the mimic system, a seemingly very reasonable requirement. The solution of equation Q has been shown 
to yield accurate structure for both ionic solutions Il22ll29ll and molecular water ll2lll in nonuniform systems, and a simple linear 
response method for solving the above equation has been derived ifisll . leading to fast and computationally efficient solutions of 
the LMF equation. 

Site-site pair coiTelations in bulk fluids may be simply related to those arising from fixing a given site at the origin, thus 
allowing us to describe structure in uniform fluids from the nonuniform perspective of LMF theory in equation (|4]i. As such, 
for bulk molecular fluids with spherically symmetric site-site interactions as considered here, we would expect that the general 
LMF equation ^ in this case could then be written as 

VK|,(r) = ^^.oM + i j dr'p?s,e.(r'|OH(|r-r'|), (5) 

with Tj being the site fixed at the origin as indicated by the conditional notation \rj on the left side of equation (|5]). In analogy 
with equation the first term is the short-ranged potential due to the only fixed charge in the system, the charge from site -q at 
the origin. The charge density p|j tot(^'l^) ioXsA charge density in the nonuniform mimic system with |0 again indicating 
the fixed site. In the case of these small site-site molecules, this total charge density may be decomposed as 

• the intramolecular charge density arranged around the molecular site 77 fixed at the origin including the site r], denoted 

^'fl.,J\/|^(^|0)' and 

• the intermolecular charge density from other unconstrained mobile molecules induced by i], denoted p|j(r |0). 

While only 77 contributes to the Vo in this equation, the intramolecular sites attached to the site i] fixed at the origin contribute 
directly to the total charge density tot ^^'^ ^1^° implicitly but strongly impact the form of the intermolecular charge density 
based on their inclusion in the simulation of the mimic system. 

While the discussion above should make the form of equation (|5]) quite plausible, two further approximations are needed in 
addition to the three stated in Appendix|A]to carefully separate effects of intra- and intermolecular charges in this equation and 
to assess its accuracy for site-site molecular models. Again, we employ the exact relationship between the pair distribution 
functions in a uniform fluid and the conditional singlet density profile due to a site fixed at the origin. 

The YBG hierarchy for site-site molecular systems with one site of one molecule fixed at the origin is derived in AppendixlB] 
The derivation is quite interesting technically, since we use an external field to localize only a particular molecular site at the 
origin rather than to represent an entire fixed molecule, as is usually done. Moreover, we derive first the YBG hierarchy for 
correlation functions between specific molecules rather than the usual generic correlation functions used in standard treatments. 
These features allows us to more easily disentangle contributions from intra- and intermolecular correlation functions. Using this 
new YBG hierarchy, the derivation of LMF theory for a uniform fluid of site-site molecules then follows the traditional route, 
while requiring two new but very plausible approximations related to intramolecular correlations, as shown in AppendixICl This 
provides a rigorous derivation of equation (|5]l. 

Towards the goal of understanding the accuracy of SCA truncations in uniform fluids, we rewrite equation (|5]l in a way that 
focuses on the long-ranged contributions to Vr, which we call Vri. 

V«il,(r)^V^I,(r)-^t;o(r) = i| drV?j,tot(^'|OM (|r - r'|) . (6) 
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If Viii\rj ~ then simulating with spherical truncations alone as in the SCA will give very accurate results. 

This Vflii^ defined in equation (|6]l is seen to be the Gaussian-smoothed electrostatic potential arising from the total charge 
density in the fluid induced by the charge from fixed site 77. This total charge density includes the single molecule charge 
distribution j\/|^('"|0) well as contributions from other fully mobile molecules. As discussed in Refs. ITtIi . ||2i1i . and 12911 . 
the restructured electrostatic potential V^i (r) induced by a general fixed charge distribution pj^^^ (r') satisfies the single Coulomb 
LMF equation given by the convolution of vi (r) with tot('^')' including contributions from both fixed and mobile charges, so 
equation ^ has exactly the form that would be expected. 



IV. SUCCESS OF SCA EXPLAINED 

We specifically explore the meaning and consequences of the LMF equation for SPC/E water, both because it has fewer sites 
than acetonitrile and also because it has a fixed geometry, thereby allowing for analytical determination of A/|r)(^l^) without 
simulation and independent of perturbations from other mobile molecules. For either hydrogen site, 

(^R,M\H(r\^) = lHHr) + qo^—2 2 ' (7) 

and for the oxygen site, 

oImioW = loS (r) + 2qJ-^^—P^, (8) 

where the charge densities have been spherically averaged about the site fixed at the origin. Separating out the contribution of 
these intramolecular charge densities, the total Vj^i|^ in equation (|6]l may be decomposed into intramolecular and intermolecular 
contributions as 

= 7 / ^^'4A/h('^'|0H'i(l^-^'l) + 7 / dr'p%{r'\0)v,{\r~r'\), (9) 

where the first term corresponds to the long-ranged interactions due to sites within the molecule with site rj at the origin. 

The success of SCA shown in Section|II]suggests that V^^ii,, « is a well-founded approximation. Before utilizing the ana- 
lytical charge densities above, we first explore an alternate formulation of LMF theory for site-site molecules which might seem 
initially fruitful. Theoretical treatments of molecular models often involve fixing a given molecular orientation and considering 
the fluid response to this configuration. Based on the splitting of 1/r, the majority of the strong electrostatic potential energy 
and force will be included in the t'o(r) used as a pair potential in SCA simulations. However, for any one orientation of a water 
molecule, the combined forces due to the vi (r) on other oxygen and hydrogen sites will not be negligible, even though they are 
slowly-varying on the scale of a for each individual site. 

As one example of this, the long-ranged electrostatic potential arising from a fixed orientation of a water molecule with O 
at the origin and the 1/r interactions replaced by vi{r) with cr = 4.5 A is shown in Fig.|6] Based on this single snapshot of 
the vi (r) contributions due to intramolecular sites, neglect of these long-ranged forces in the SCA would seem an ill-conceived 
approximation, and we might suppose that a Vr depending on both intermolecular distance and relative molecular orientations 
would be required. However, looking at an individual orientation of the water molecule for long-ranged interactions fixes all 
three intramolecular charges and would be expected to generate a very different and larger density response than the single fixed 
molecular charge at the origin needed to determine radially-symmetric site-site correlation functions, as the LMF equation ^ 
and Appendices iBl and ICl show. 

The first term in equation (|9]l may be determined analytically for SPC/E water, and this is the first crucial step in understanding 
why the full Vj^i\n will be small to good approximation in uniform systems. As shown in Fig.|2l the spherically symmetric long- 
ranged potential from the first term, which we shall term V^ji intral?? (^) is indeed much more slowly-varying than the orientation- 
dependent potential shown in Fig.|6l In this figure, we compare V^i intra|i7(^) to both the VQ{r) due to either O or H fixed at the 
origin as well as the vi{r) due simply to that site rj. 

From these plots, we see that V/fi,intra|r) is substantially smaller in magnitude than either the electrostatic potential arising 
from a specific water molecular orientation or the potential due simply to the charge at the site we have fixed at the origin. 
Therefore, we infer that the spherical truncations prescribed by LMF theory and the associated mean-field averaging of long- 
ranged interactions will actually be even more effective in bulk molecular simulations than in a corresponding simulation of 
an ionic system with charges not bound into neutral molecules. Again we emphasize that this spherical averaging is not an 
unfounded approximation, but that it arises rigorously from the statistical mechanics of molecular models interacting via site- 
site potentials. 

The total electrostatic potential arising from the spherically-averaged intramolecular charge density will be exactly zero for all 
r > Ion if oxygen is fixed at the origin or for all r > Zhh if hydrogen is fixed at the origin. Thus it might seem counterintuitive 
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FIG. 6. The long-ranged potential in the z=0 plane resulting from a fixed orientation of a water molecule with ro ~ (0, 0, 0), rni = (1, 0, 0), 
and ru2 = (—0.334, 0.943, 0) and a chosen as 4.5 A. The potential is displayed in units of ksT/eo in order to aid in gauging the magnitude 
of this potential relative to thermal fluctuations. The chosen orientation of the water molecule is shown with solid black lines and points. 



that the corresponding V^i intia|r) is small but non-vanishing beyond this distance. However, the distinct treatments of the short- 
ranged and long-ranged parts of 1/r using Gaussian convolutions in LMF theory require just such a nonzero potential. All the 
short-ranged parts of 1 /r are treated explicitly via vo{r) positioned around each site in the water molecule in order to represent 
local correlations; the capture of these local correlations in the SCA simulation is crucial. In tandem, only the long-ranged 
components vi (r) are spherically averaged about the fixed site 77 in LMF theory, leading to a non-zero but slowly-varying and 
small magnitude potential V/ji^intralijl^") outside the total potential cutoff. For the correlations between molecules, the need for 
non-zero short-ranged site-site vq terms seems quite natural; the need for similar short-ranged terms also holds for the far-less 
intuitively-obvious splitting of the (exactly zero) electrostatic potential between two charged plates ll29tl . 

As demonstrated in Fig.|2l V^^i i„ti.a|r;('') is quite small and slowly-varying for SPC/E water However, while this may make 
the approximation that the total Vi^i|^(r') « in equation (|9]l plausible, it certainly does not guarantee it. Therefore, we also 




r (A) r (A) 

(a)Large Scale (b)Focus on V^i^intralr, 

FIG. 7. Comparison of Vfli.intra|i7('') to relevant potentials due solely to the site 77 fixed at the origin, whether it be oxygen (red) or hydrogen 
(blue). This electrostatic potential due to the whole neutral molecule is substantially smaller than both the short-ranged and long-ranged 
components of 1 /r due to the individual site fixed at the origin. 
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FIG. 8. Estimation of V^jii^ based on charge densities from the simulations conducted using Ewald summation. 



estimate V^i|j,(r) by directly inserting the charge density resulting from the Ewald simulation into the LMF equation. The 
sole care we take is in enforcing overall charge neutrality at the cutoff radius for the potential, as this is also the furthest radius 
at which g{r) is calculated. As seen in Fig. [8] these potentials are also quite slowly-varying, lending strong credence to the 
approximation Vjji\ri ~ Oin determining structure. 

This approach for determining Vjin^{r) requires a full Ewald simulation, contrary to the general philosophy of LMF theory, 
which seeks to use simulations only in the mimic system. Thus strictly speaking we should self-consistently solve for Viji|^ 
based on charge densities from the short-ranged mimic system using the linear-response treatment developed in Ref. ifisll . But 
previous work has shown that the full LMF theory gives excellent agreement with the results of Ewald simulations for water 
even in nonuniform environments, so this Ewald determination should be very accurate. Furthermore, care should be exercised 
with the fc = component of any charge density used in the LMF equation lfl6l[30ll . However, charge densities obtained via 
Ewald summation exhibit exponential screening and strictly enforce overall neutrality, thus easing the need for great caution in 
the treatment of the fc component. 

This simple estimate based on the Ewald charge density certainly suffices to demonstrate that Vi^i|^(r) is small and slowly- 
varying in this case, and provides strong justification for the accuracy of the SCA. In general we expect that quick estimates 
of Vr using Ewald charge densities when such simulations are computationally practical will be very useful in obtaining an 
accurate initial estimate of the final self-consistent Vr, and one that will be almost certainly in the linear regime where the 
method of Ref. fisl] will be especially easy to use. 

However, for these bulk fluids an accurate Vri is neither necessary for determining the structure to the accuracy shown here 
nor for determining the thermodynamics of the fluid as shown in Ref. ||23ll . Provided that a sufficiently large a is chosen, simple 
spherical truncations in simulations coupled with thermodynamic perturbation theory yield accurate structure, energies, and 
pressures. In the case of SPC/E water, structure might indicate that any a > 3.0 A is sufficiently large, but thermodynamics via 
perturbation theory showed that a > 4.0 A is required 123[| . 

In general, the choice of a sufficiently large a is crucial for the accuracy of LMF theory. For the acetonitrile system at the 
higher temperature and lower density, inclusion of a self-consistent Vri with cr = 4.5 A gives a poor description of the structure 
of the acetonitrile system. However, since our simple scaling analysis suggests that tT,ni„ ~ 6.5 A, we do not expect LMF 
theory with the smaller a to be able to correct this structure. For the acetonitrile systems at low and high temperature, just as 
for the water system at ambient temperatures, a sufficiently large cr yielded accurate results simply via SCA. Furthermore, the 
acetonitrile results already demonstrate that tr does not need to be on the scale of the entire molecule but rather on the scale 
of nearest neighbor correlations, as is expected from derivations of LMF theory. In Appendix iDl we discuss LMF theory for 
CHARMM-like molecules in order to better state the necessary conditions for choice of a in much larger molecules. 



V. CONCLUSIONS 

In this paper, we have demonstrated the accurate results possible using spherical truncations of 1/r interactions in simulations 
of uniform fluids. We show that these spherical truncations yield not only highly accurate pair correlation functions but also 
highly accurate dipole-dipole coiTelation functions. This good performance in bulk simulations of pair correlation functions was 
known; however, a solid theoretical justification for the use of such spherical truncations in molecular systems has been lacking. 
In this paper, we present just such a theoretical backing - local molecular field theory. The derivations relevant to LMF theory 
for a variety of site-site molecular models are presented in appendices and the main paper focuses on understanding the accuracy 
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of these spherical truncations both phenomenologically and quantitatively using LMF theory. LMF theory provides a general 
conceptual framework that helps us understand why spherical truncations generally work so well in uniform systems and also 
provides the essential corrections needed in most nonuniform environments. 
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Appendix A: Yvon-Born-Green (YBG) Equation and Local Molecular Field (LMF) Theory Derived for Small Site-Site Molecules in 

an External Field 

While this paper primarily deals with a uniform site-site molecular fluid, the derivations of both the YBG hierarchy as well 
as the LMF equation are simpler for a general nonuniform system. We present this derivation here using a straightforward 
method that also introduces the basic site-site notation and ideas that we will then generalize and apply to the uniform fluid 
in Appendices IB] and ICl where careful attention is paid to the distinction between intra- and intermolecular correlations. We 
should note that MuUinax and Noid 13 lil have developed a basis expansion method that can be used to derive a generalized YBG 
equation for a variety of molecular systems. 

Previous work 13211 for site-site YBG equations begins the derivation by writing the singlet density for a molecular site in 
terms of the singlet density for the entire molecule with a fixed orientation, taking appropriate gradients on either side, and only 
then reducing to a site-site representation. Using the general formalism developed by Chandler and Pratt ||33ll for the partition 
functions and density distribution functions of mixtures of site-site molecular models, we may follow a similar path to the 
derivation of a general site-site YBG equation. The formalism originally was developed to also account for the possibility of 
chemical reactions, and since this is not a concern in the inherently classical systems we study, a few alterations will be made to 
simplify notation, with no impact on the meaning of the equations. 

The partition function for a mixture of molecular species M with total sites um on each molecule labeled by Greek characters 
such as ^ is given below with the position of the ^ site on the i'^ molecule of type M given as r|^^ and the positions of all um 
sites on the molecule of type AI given as RiM- 

Q{{M}) = i\{N,A4r n {^mT^'\ I {\{d^-.^ (Al) 

where the total potential energy U is defined as 

Nm Nm riM 



2^ = E E ^M(RiM) + E E E M^Z) 

M i=l M i=l ^=1 

^ Nm Nmi riM n^,, 

+ 2 ^ ^ ^ ^ ^ <5a/M'(5^j ) E E ^iMaM' 
M M' i=l j = l 5=1 Q=l 



(A2) 



Here I'm is the symmetry number of the molecule. For example, for H2O, i/ ~ 2 for 2 equivalent orientations, and for CH4, 
1/ = 12 for 12 different equivalent orientations - 3 equivalent rotations for each of 4 different C-H bonds fixed in position. With 
symmetry numbers included, each "equivalent" atom may be correctly viewed as a unique site. Thus H2O has 3 sites and CH4 
would have 5 sites, h!'^'! is the thermal de Broglie wavelength for the atom ^ on molecule M. The factor of {\ — 5MM'5ij) ensures 
that the general pair interactions u^moM' , often taken as a sum of Coulomb and LJ interactions in CHARMM-like models, arise 
only for sites on different molecules. We will consider modifications necessary to apply this reasoning to a true CHARMM 
model for larger molecules in Appendix iP] 

We now write the single-site density distribution function using the notation dR to represent all molecular coordinates Ra/, 
and "division" by dr'^^^j to indicate integration over all particle positions except the f site on the molecule of type M. Thus, 
we have 

.(1) (r.\-^ f .-PU ( ^ 

dv 



P^Mr) - ^ ) , (A3) 
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with Zn the configurational partition function and normalization constant given by integration over all R. Here, r has replaced 
'^iM ^- Now taking the gradient with respect to r and using the equivalence of all molecules of a given type, 

dR \ 



dr 



(0 

IM, 



[V0M,c(r)] 



Nm 



-/3U 



dR 



Nm{Nm - 1) 
Z 

NmNm' 



r — r- 



(a) 
2M 



M'^M Q=l 



lA/ / 

dR 



(A4) 



M , 



We may simplify this site-site molecular YBG equation in terms of an intramolecular density distribution function, pa/(Rm), 
and a two-point intermolecular site-site density distribution function, p^lj„j^,j, (r, r')' specifically defined to exclude intramolec- 
ular site-site correlations. Here we set Rim = Rm, ~ ^2m' ^ ^' ^- 



eA/i-tV-A/j — / 



/ dR 



VdR 



lA/ 



„(2) - Nm{Nm' -Smm') f pu 



dR 



(A5) 
(A6) 



Substituting these definitions into equation ( IA4b . we find 



dR 



M 



dr 



ii) 

M 



[V0A/,?(r)]pg,(r) + ^^ / drVgU,,(r,r')Vii5A/.A/'(|r-r'| 



(A7) 



A/' a=l 



The sole difference between this equation and the YBG equation for atomic mixtures is the term involving the gradient of the 
bonding energy and the intramolecular density distribution function. In order to put this exact YBG equation in a standard form 



from which the LMF equation is derived, we divide each side by p^^f (r), yielding 



fcijTV (in pjy^ (r)) = J [Va;A/(RM)] tok(R-Mlr) 



dR 



M 



dr. 



(0 

M 



"m' p 

+ V(/)A/c(r) / dr'paA/,|^A/(r'|r)Vit^A/aA/'(|r - r'l). 

M' a=l-' 



This division generates conditional densities on the right side of equation ( IA8I ). Thus 

.(2) 



(A8) 



(A9) 



is an intermolecular conditional density, proportional to the probability of finding site a on a molecule of type M' at position r' 
given that site ^ on a molecule of type M is located at position r, and similarly I^IR-mIt) is the intramolecular conditional 
density of a molecular orientation Rm given that site ^ is located at position r. 

We now derive the LMF equation. We first consider a general separation of the intermolecular interactions into short- and 
long-ranged parts 



l-aM^M'{r) = Uo,qA/C A/' (r) + Ui^aM^AI' (r) , 



(AlO) 



where ui is slowly varying over the range of strong nearest-neighbor interactions. We seek a mimic system which is composed 
of molecules with only short-ranged intermolecular interactions UQ^aM^M' (r) along with effective single-particle potentials 
4>R,iM{Y), chosen in principle so that the induced singlet densities in the full and mimic systems are equal: 



(1) 



(All) 
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All intramolecular and bonding potentials will be assumed to be the same in the mimic and full systems. 

Following the standard path to the LMF derivation, we take the exact difference between the YBG equation for the full system 
and the YBG equation for a mimic system, assuming the equality of the singlet density profiles. After rearrangement we find 

"a/' p 

Viti^^MaA/' (|r — r'l) 

A/' 4=1 

+ I {£'A/|4(R-M|r; [0]) - Qii^M\i{'^m\Y] [Vwa/(Rm)] ( | 

+ / '^i"' {PaA/'|af(i"V; ['/']) - Pfl,aA/'|a/(r'|r; Vuo,4A/aA/'(|r - r'l) 

A/' Q = l 
"aj' !■ 

+ rfr'{p„M'|4A/(r'|r;[0]) -PaA/'(r';[0flJ)} Vui,4A/aA/' (|r-r'|). (A12) 

A/' a = l 

The above equation is exact but not particularly useful as it stands because of the appearance of complicated conditional 
densities on the right hand side. In order to yield the LMF equation, we must make three connected and very reasonable 
approximations for the integrands of the last three terms based on our chosen forms for ito and ui. 

• Approximation 1: The densities of specific molecular orientations will be well approximated by the mimic system such 
that 

fi-AfislRivilr; [0]) ~ £»i^^A/|^(RM|r; (A13) 

allowing neglect of the integrand involving these functions. For small molecules, this seems like an eminently reasonable 
approximation, since the prevalence of various relative intramolecular orientations in both systems will be dominated by 
the identical short-ranged interactions and the overall molecular orientation should be quite well approximated given local 
short-ranged interactions and the long-ranged orientational corrections due to Vj?. 

• Approximation 2: The product {pQA/'|eA/(r'|i"; - PR.aM'\^M{r'\r; [(pn])} Vuo,4a/qA/' (|r - r'|) can be neglected. 
This term probes the difference between the conditional singlet densities for the full and mimic systems via convolution 
with Vuo{r). The integrand will be quickly forced to zero at larger |r — r'| by the vanishing gradient of the short-ranged 
uo{r). The integrand will also be negligible at small |r — r'| since both the full and mimic systems have the same strong 
short-ranged core forces with an appropriately-chosen Uo{r), so the density difference inside the curly brackets should 
then be very small. 

• Approximation 3: The final product {paA/'|4A/(r'|r; [(/>]) — PaA/'(i"'; [0_r])} Vui^^MqA^' (|r — r'|) can also be ne- 
glected. This is due to the fact that difference between the conditional singlet density and the singlet density of the full 
system will be most substantial for exactly the small distances where ui is slowly varying and Vui (|r — r'|) will be 
small. At large separations the conditional singlet density reduces to the usual singlet density except in special cases like 
near the critical point, so this term can again be neglected. 

Approximation 1 is the sole new addition as Approximations 2 and 3 are identical to those required for single site mixtures as 
detailed in Ref. ITtIi . However, when these reasonable approximations are employed and LMF theory is applied only to the 
charge-charge interactions of molecular models so that charge densities can be introduced as in IitIi . we can exactly integrate 
the remaining terms in equation ( IA12b and find the desired LMF equation for site-site molecular models: 

(pR.^Alir) = 0rie,CA/(r) + 9CA/Vfl (r) 

Vi?(r) = V(r) + iy"drV|j(r')i'i(|r-r'|). (A14) 

Here (/)„e,4A/(r) contains all the non-Coulombic parts of the external field and V(r) is the electrostatic potential from the fixed 
charge distribution as explained in detail in [ITfl. Each molecular site now moves in a renormalized electrostatic potential 
due to an average charge density p|j(r) that is partially contributed to by it and its bound molecular sites. This might seem 
to be a cause for concern, since implementations of Ewald summation do remove the effect of both the charge itself and these 
bound charges [[34}]. However, we argue that this is reasonable since LMF theory convolutes the average charge density, not the 
instantaneous charge density, with the slowly-varying long -ranged ui(r). 

The equation above is identical to the mixture LMF equation as related in previous derivations. However, the preceding 
derivation for small site-site molecules helps us to understand that the use of the mixture LMF equation for site-site molecules 
still is grounded in the YBG equation with solid statistical mechanical approximations. It also sets the stage for the notationally 
more complex derivations for bulk site-site molecules given below. 
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Appendix B: Derivation of YBG Equation Appropriate for Uniform Small Site-Site Molecules 



Now we derive the YBG equation for pair distribution functions in a uniform system of small site-site molecules. We first 
consider a system of only one molecular type in order to focus on the new features needed to easily separate out contributions 
from intra- and intermolecular interactions. It is straightforward to generalize these results to a mixture of molecular types as 
indicated at the end of the appendix, and this method can also provide an alternate derivation of equation ( IA8I ) as well. 

Our broad strategy in deriving the YBG equation for site-site pair distribution functions in uniform fluids uses the equivalent 
functional forms of the pair density function and the conditional singlet density. The conditional singlet density may be physi- 
cally interpreted as the density that would arise if a single site were fixed at the origin. In the following derivation, we apply a 
special external potential in the Hamiltonian which yields exactly this situation. Note that this is different than the standard use 
of an external potential to represent an entire molecule with a given orientation fixed at the origin. Due to the intramolecular 
correlations, several new terms arise in the YBG hierarchy. 

In a classical system even identical molecules or sites can be treated as distinguishable. It will prove useful to generalize the 
external fields 0'") appearing in the Hamiltonian in AppendixlAlbv assuming that the system interacts with a set of external fields 

— that in principle can differ for each site a of each molecule i. The total potential energy of the nonuniform 

molecular system with this very general set of external fields can then be written as: 

N n N ^ N N 



i—l a—1 i—1 2—1 j—1 a—1 ^—1 

We first consider molecule-specific distribution functions like 

Pf^(rP;[0]) = ^/e-^"'^)^, (B2) 

the probability density for finding site ^ of particular molecule 1 at r and will later consider the usual generic distribution 
functions like that given in equation ( IA3b . which account for the equivalence of molecules of the same type. By taking the 
gradient of equation ( IB2I ) we immediately derive the first equation of the specific YBG hierarchy: 



Pm(Ri;[0])V ,oWA/(Ri) 



N 



J V It n 

+ / P^l\r^^\rf-Myrioua\r['^-rf\)drf. (B3) 

This YBG equation is identical to that derived in Appendix |A] with the important difference that it does not appeal to the 
indistinguishability of molecules of the same type. This is crucial because the external field we will apply explicitly fixes one 
site on a given molecule at the origin. Here Pm (Ri', [</>]) in the second term on the right denotes the n-site intramolecular 
distribution function, defined as in equation ( IB2b but with integration over Ri excluded. The integration in the second term is 
over all Ri with site ^ fixed at r^''''. Similarly the definition of P^'^ {r^i \ r^"''; [0]) excludes integration over r^^' and r^""* and 
involves sites on different molecules 1 and j. 

We want to determine intermolecular site-site pair distribution functions in the uniform system with = 0: P^"^ {r'"P , r^"'' ; [(j> = 

0]) = P(a\\''^''P ~ ''i"' I)- Even with a general anisotropic ujm, these can depend only on the radial distance between sites ^ and 
a on different molecules 1 and j because of translation invariance and the spherical symmetry of the intermolecular potential 
Mq.{ in equation (IBlb . 

We gain information about these uniform system functions by considering another special case of equation ( IB3b where only a 
single field (j)^2^ (rj''' ) involving a given site i] on a particular molecule 2 is nonzero. This field has a special form that confines 
this site to a very small spherical region centered about the origin 0. Thus 4>^2^\r^2^) = oo for |r^''''| > e and is zero otherwise 
and we are interested in the limit e — s> 0+. All other are zero. In order to aid in visualization of the various sites and 
molecules, the basic inter-relation of site indices used in this appendix are shown in Fig.|9] 

Note that the nonzero field 02^'' {r'^^ ) only appears implicitly in equation ( IB 3b through its effect on the distribution functions 
and that this field fixes only the single site rj of molecule 2 at the origin, and not the orientation of the entire molecule. In the 

limit e 0+, P^^\r^^'^; [(f>'^'']) in equation ( IB3b reduces to a conditional singlet density with site rj of molecule 2 fixed at the 
origin. Taking account spherical symmetry we write this as 

pf)(rP;[4'')])=P«(r|0), (B4) 
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where we set r ~ and note that P^^J^ can depend only on the magnitude r of r. The bar before the subscript i] and 
the argument on the right side indicates a conditional density with site 77 in constrained molecule 2 fixed at the origin. By 
translational invariance the specific pair distribution function P^^^ (r) in the uniform system equals times the corresponding 
specific conditional singlet density in equation (IB4t . and we will use this equality later to determine uniform system pair 
distribution functions. 

The nonuniform pair distribution functions in equation ( IB3I ) can be similarly rewritten in this special case. In particular, the 
pair distribution function involving another site 7 ^ 77 on constrained molecule 2 can be written as 

p(^)(r(«\4^);[4'"]) = ^S,(r,r'|0) (B5) 

where we set rj^' = r'. In this and the following appendix we will generally use a single prime to denote coordinates on 

(2) 

the constrained molecule. P^^^^j is strongly affected by the fixed site and the short-ranged intramolecular interaction ujm in 

equation ( IB 11 1 and vanishes for large |r'|. This is even more true for the distribution function P^^^'^ (r^^'' , r2'''' ; [02 '''])' which has 
the hmiting form as e — 5- 0+ 

P(f (ri^^r^"); [4")]) = P^ (r|0)<5(r('') - 0). (B6) 

Both these distribution functions are very different from those involving any site a on an unconstrained third molecule, which 
takes the form 

pg'(r(«\r("M4")])=p(^j,,(r,r"|0) (B7) 

where we set rg"'' = r", and generally use double primes to denote coordinates of unconstrained molecules. In contrast to 
equation ( IB5I ), this does not vanish for large |r"|, where it reduces to a product of conditional single particle functions for large 
|r-r"|. SeeFig.m 

We also define an induced single particle interaction on site ^ associated with the pair potential from fixed site 77 at 0; 

05l^(r) =U5^(r) (B8) 
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and rewrite equation ( IB3l l using the new notation in this special case. Separating terms involving constrained molecule 2 from 
those that involve other unconstrained molecules, we get 



/dR. 
PM|r,(R|0)VrWAf(R) ^ 

n „ 

+ P^'iir,r'\0)VrU^,{\r--r'\)dr' 

a = l 

Using Eqs. ( IB6I ) and ( |B8t . the first term on the right side of equation (IB9I ) is generated by the a = 77 and j = 2 term in 
equation ( IB3I ). where the pair potential from the fixed site 77 acts like an effective external field on site ^. We have used the 
equivalence of all molecules except 1 and 2 in the last term in equation i 



To get to the final form useful for LMF theory we divide by P^^^{r\0) and introduce the usual generic distribution functions. 
Thus the distribution function for finding site ^ of any other molecule at r is 

p[]l{r\0)^{N-l)P^l^^{r\0) (BIO) 

Similarly the generic distribution function involving three distinct molecules in the last line of equation ( IB9I ) is 

pfj|^(r,r"|0) ^{N- l)iN - 2)P^l]^ir,r"\0) (Bll) 

Division by p^|^(r|0) will yield a density conditioned by ^ as well, defined by 

pW^(r"|0,r) EE p[^j|,^(r,r"|0)/pW(r|0) (B12) 

The remaining distribution functions in equation ( |B9t involve sites on only two molecules and have very different forms 
strongly influenced by the intramolecular interaction ujm- We again use the symbol g to emphasize this point and define generic 
functions 

eM|„(R|0) ^{N- l)PM|r,(R|0) (B13) 

and 

4;'|,/r, r'|0) ^{N- l)P^X^r, r'|0). (B14) 

Densities conditioned on ^ as well are similarly defined as in equation ( IB12I ). 

Using this notation in equation iB9i we arrive at the desired final form for the first equation of the site-site molecular YBG 
hierarchy, with site ?/ of a particular molecule fixed at the origin: 



fcBrVrlnpJ|,5(r|0) = Vr05|^(r) 

gM|„c(R|0,r)VrCJA/(R)^ 



n « 

^ J 4;i^(r'|0,r)Vr^.^^(|r-r'|)dr' 

n „ 

E J pi|!,5(r"|0,r)Vr7.c4|r-r"|)rfr". (B15) 



Note that this YBG equation is nearly identical to the equation derived in Appendix lAl An important additional contribution 
arises due the correlations between the site ^ and the various sites 7 present on the same molecule as the site ?/ fixed at the origin. 
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It is straightforward to extend this approach to a general mixture of molecular species. With obvious generalizations of 
notation we find for site ^ of species AI with site 77 of a different molecule of a possibly different species M' fixed at 0: 

+ 4M'hAf'«Af('''|0'i")^r"a/7M'(|r-r'|)dr' 

+ EE / pi'i,.|,M'eA/(r"|0,r)Vr^/5MoA^"(|r-r"|)rfr". (B16) 

Af" a=l 

Appendix C: Derivation of LMF Equation Appropriate for Uniform Small Site-Site Molecules 

We now derive the LMF equations appropriate for a uniform mixture of site-site molecules, using the exact YBG equa- 
tion ( IB16b . The basic strategy follows that for the molecular system considered in Appendix|A] We again consider a general 
separation of the intermolecular interactions into short- and long-ranged parts, as in equation jAlOl ) such that the mimic system 
will have only short-ranged intermolecular interactions along with effective single-particle interactions (f>R^^M\7]M' {">') associated 
with the fixed site at the origin. These effective interactions are again chosen in principle so that the induced densities in the full 
and mimic systems are equal: 

pSm|.a/'(^|0)=4m|,m'(H0). (CI) 

All intramolecular and bonding potentials will be assumed to be the same in the mimic and full systems. In Appendix iDl we 
generalize to instances of larger molecules where long-ranged interactions might exist between sites on the same molecule. 

Following the standard path to LMF derivation, we take the exact difference between the YBG equation for the full system 
and the YBG equation for a mimic system in a restructured field for which equation ( ICll ) holds. Since we already must include 
subscripts for the fixed site and two other sites, for simplicity of notation we will first consider a single component site-site 
molecular system. 

Thus, using equations ( IClb and ( IB15I ) we have exactly 

Vr[0fl,5|^(r-) = J |eM|r,5(R|0,r)-efl,M|„?(R|0,r)|vrCJA/(R)^ 

; J |p(y,,^(rnO,r)-pW„|^^(r''|0,r)|vr«o.?o(|r-r''|)dr'' 

+ E / |pS,,(r"|0,r)-pi^;(r"|0)|vr^i,5„(|r-r"|)dr" 

71 p 

n „ 

+ E / P?o|.(^"|0)Vr^i,^„(|r - r"|) dr" (C2) 

As a consequence of our judicious choice of ito(r) and ui(r), all the integrals involving terms with large curly brackets 
vanish to a good approximation. The first, third, and fifth terms with curly brackets all may be neglected by Approximations 
1-3 as detailed in Appendix lAl We now must employ two related approximations leading to cancellation of the intramolecular 
correlations functions. 



n 

E 

a=l • 
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• Approximation 4 The product |p^"'||^^(r'|0, r) — p^''^|^^^(r'|0, r)| VrMo,C7(|r — r'|) will be approximately zero. The 

logic here is virtually identical to that of Approximation 2. Given rigid or even flexible bonds between intramolecular 
sites 7 and ry we expect the matchup between densities in the short-ranged system and the full system to be even better at 
short distances, leading to an even stronger cancellation. 

• Approximation 5 The product |£i^^^]^^(r'|0, r) — gi^|^(r'|0)| VrUi_^^(|r — r'|) will also be approximately zero, for rea- 
sons similar to Approximation 3. In fact, the intramolecular conditional density profiles should be less sensitive to the 
presence of a site on another molecule for many configurations. At small separations, the cancellation due to the slowly- 
varying nature of ui{r) will still hold. 

Thus we see that while more intramolecular terms must cancel in the derivation of the LMF equation, exactly the same fine 
of logic is followed as in Appendix |A] Using (IClb and setting the first five integral terms to zero as justified in the previous 
discussion, we arrive at the site-site LMF equations for each combination of fixed site t] and mobile site ^: 

<I^RM^)-khi^)^Y. J ^'?-y|,('^'|0K.C7(|r-r'|)dr' 

+ E / p£|,('''|0)ui,,„(|r-r"|)dr" + C. (C3) 

a=l-' 

In the above equation, there are terms due to intramolecular sites as well as sites on other molecules. The portion due to 
intramolecular sites does not imply an action of 0^ on intramolecular sites but rather includes the effect of these intramolecular 
sites on sites of other molecules. This set of equations for each choice of ^ and rj has the simplest form possible with a general 
separation of the pair interactions into short- and long-ranged parts. 

However, as discussed in detail in iITtIi . LMF theory takes a particularly simple and powerful form when it is applied only 
to Coulomb interactions and all charges are separated using the same cr, as we do in this paper Charge densities rather than 
individual molecular site densities can then be naturally introduced, as shown below. Furthermore, these new equations based 
on charge densities are not only simpler, but likely lead to an even stronger overall cancellation of terms than argued for each of 
the individual terms previously. 

We may write the long-ranged Coulomb part of the specific intermolecular pair interactions as 

ui,aMiM'(r) = vi{r) (C4) 

and as before the short-ranged core interactions will be defined as UQ^aM^M'ir) — UaM^M' (r) — ui^aM^M'{r) and will en- 
compass all LJ-like interactions as well as the usual Coulomb core Vf){r) terms. In particular, using equation ( |B8t , the induced 
interaction from the fixed site can be written as 

'I'ilvi'^) = <^ne,^|r,(^) + +Vi{r)], (C5) 

where 4>ne,^\r]{'r) contains all non-electrostatic (usually LJ) pair interactions between sites ^ and rj. 

The relevant spherically symmetric charge distribution arising from the constrained molecule with site i] fixed at the origin is 
given by 

n 

SR,M\v^r\0) ^ q,S{r - 0) + ^ g^g^ |^(r|0). (C6) 

For rigid molecules like SPC/E water, a/|,,(''|0) can be determined in advance and expressed solely in terms of sums of 
(5-functions as discussed in equations O and 

In general, the total induced equilibrium charge density for a site rj fixed at the origin is then 

n 

Pho^O) ^ g%,,^^^ir\0) + ^ qj^l^,.^{r\0), (C7) 

Q = l 

where the second term is the contribution to the charge density from the other unconstrained mobile molecules. 
Using (IC4b . we see equation ( IC3b can now be written in the compact form 

(t^RMvi^) = (t>ne,^\ri{r) + 9? Vi?,|,, (r) , (C8) 
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where Vfl|,,(r) is the restructured electrostatic potential induced by the fixed site rj and the other associated intramolecular sites. 
This satisfies the site-site Coulomb LMF equation 

VRi.ir) = ^vo{r) + \j dv' pl,^,{r'\G)v, (|r - r'|) . (C9) 

We may also define a restructured potential Vfli|,, containing only the long -ranged components of the potentials as 

Vfli|,(r) = Vfl|,(r) - '^vo{r) = -J dv' p% ,^,{r'\0)v, (|r - r'|) . (CIO) 

This is the restructured Gaussian-smoothed electrostatic potential induced by the fixed charge from the site 77 at the origin, where 
Pr tot ^oi^\ equilibrium charge density due to the fixed charge, the intramolecular charge density, as well as the fully mobile 
charges on other molecules. 



Appendix D: Treatment for Non-Uniform and Uniform Larger CHARMM-Iike Molecules 



While the use of Approximation 1 in previous derivations of LMF theory might suggest that our findings are invalid for 
larger molecules defined by CHARMM- or AMBER-like potentials, this is not the case. Far less restrictive approximations than 
the equivalence of the whole molecule density functions may be derived by the usual postulates of a bonding potential form 
where only sites separated by 1, 2, or 3 consecutive bonds in a molecule may experience a special local bonding interaction. 
Pairs of intramolecular sites with larger separations interact only through spherically symmetric (usually Coulomb and LJ) pair 
potentials. 

Rather than proceeding through logic identical to that found in the previous appendices, we instead outline the approximations 
necessary. Then we briefly describe features of the LMF equations valid for such larger molecules in a general external field and 
in a uniform fluid. We seek to emphasize that LMF theory is equally valid for large molecular models typically employed, based 
on physically reasonable approximations. 

A separate appendix dealing with YBG equations and LMF equations for larger site-site molecules is necessary because in 
most simulation potentials, such as those defined by the CHARMM ^ and AMBER parameter sets, the potential energy due to 
"intermolecular" interactions (LJ interactions and point charge interactions) is not written as distinct summations over molecules 
and their intramolecular sites. Rather the LJ and charge interaction contribution to W is a sum over all sites separated by at least 
three bonds {i.e. excluding atoms bonded or connected via angle bending). 

The expression for the partition function Q does not change, but U does. Specifically, we decompose the general lom into a 



set of bonds (^l^^j^j, bond angles w^'^^^^^^, and bond dihedrals uj'^^^g^j^j connecting appropriate sets of neighboring sites. We also 
introduce a bonding matrix Bm ct) for each species M which is 1 if sites ^ and a can be connected by two consecutive bonds 
and otherwise. With this notation we write 



.id) 



Nm riM Nm 

M i=l ^=1 M i=l a--f 

Nm 



1^ ^iM '■iM J 



EE E 

AI i—1 a — 7 — (5 
Nm 



''iMJ 

SrSr \^ /'r(") r^'^) r^-^) r^'^^ ^ 

AI i—1 a — 7 — d— 

Nm Nm' tim "m' 

M M' i=l j=l C=l a = l 



Nm Nm' um "m' 
2EEEEEE(^~ ^MM'k]BM{£.,a)) U^MaM' ( 



(Dl) 



U written in this way is virtually identical to the small site-site molecular U other than the decomposition of the bonding 
potentials and the allowance for intermolecular-like interactions between sufficiently separated sites within a single molecule. 
The first three u terms are for bond vibrations, angle vibrations, and dihedral rotations of two bonds around a connecting bond. 
Technically, these usually depend on only r, 6, and (j> respectively, but we include positions for generality and for ease in taking 
gradients in deriving the appropriate YBG and LMF equations. These sums are understood to count sets of atoms connected via 
bond, angular, or torsional potentials only once. One complication for the AMBER force field is that non-bonded interactions are 
scaled down for 1 -4 (dihedral) pairs. LJ interactions for 1 -4 pairs are divided by 2.0 and Coulomb interactions are divided by 1.2. 
We will not address this complication, but it conceivably could be included in the Bmict, 7) formalism by introducing matrix 
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elements accounting for these scalings. The new all-atom force field for CHARMM does not scale the Coulomb interactions for 
1-4 pairs. 

Based on the form of the potential, it is quite logical that Approximation 1 presented in Appendix |A] now becomes a se- 
ries of approximations related to particles connected via bonding potentials. Following the same path for LMF derivation in 
Appendix|Al we find that the weaker conditions for accuracy replacing Approximation 1 are: 

• for sites a and 7 connected via bonds, 

,(2) (rM,r(^); [0]) ^ (r("),r('^); [M) (D2) 

• for three sites a, 7, and 6 connected via a bond angle, 

,(3) (r("),r(-),rW; [0]) ^ gf (r^, r(-)r(^); [m] (D3) 

• and for sites a, 7, S, and ( involved in dihedral rotations, 

,(4) (rW,r(-),r(^),r(':);[0]) ^ gf (r^, r(-)rW , r(?); . (D4) 

These approximations are much more easily supported by mimic systems with reasonably small cr. This a may have to be on the 
order of 1-4 distances since 1-4 pairs have Coulomb interactions. In general though, we expect that LMF theory can be applied 
in standard biomolecular all-atomistic simulations with reasonable success with a a spanning only a few bond lengths rather 
than an entire biomolecular radius, as suggested by our results for acetonitrile in the main text. 

Provided that these approximations hold, we find exactly the same LMF equation for a nonuniform system. Analysis for 
the bulk uniform fluid becomes more challenging as we must in principle consider all three-particle combinations of the three 
sites - a, ^, and the fixed site 77 - where a and £, interact via their pair potential. This involves a wide range of permutations 
across different molecules, resulting in a larger number of terms that have to cancel in the derivation of the LMF equation, but 
again the final equation is essentially identical with the same underlying physical intuition. In fact, for these large molecules the 
attractiveness of treating solely charge-charge interactions via LMF theory becomes apparent. Tracking the net charge density 
profile about a site is far more manageable than tracking all possible site-site density profiles. 
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